module fouriers2d_m
	
	use MKL_DFTI

	implicit NONE
	
private

	public			:: FourierForAllSites

contains
	

pure function mymod( n, nl)
	implicit NONE
	integer(4)				:: mymod
	integer(4), intent(IN)	:: n
	integer(4), intent(IN)	:: nl
	mymod = mod(n-1,nl) + 1
end function mymod

subroutine SiteCorrelator( nl, nx, ny, sits, sitecorr)
	! Variables
	implicit NONE
	integer(4), intent(IN)				:: nl
	integer(4), intent(IN)				:: nx
	integer(4), intent(IN)				:: ny
	integer(4), intent(IN)				:: sits(nl,nl)
	complex, intent(OUT)				:: sitecorr(nl*nl)
!---------------------------------------------------------------
	integer(4)						:: ix, iy
	type(DFTI_DESCRIPTOR), pointer	:: MD_Handle
	integer(4)						:: stt
	integer(4)						:: sz(2)

	! Body of SiteCorrelator
	
	forall( ix = 1:nl, iy = 1:nl)
		sitecorr((iy-1)*nl+ix) = sits( nx, ny) * sits( mymod(nx + ix,nl), mymod(ny + iy,nl))
	end forall

	sz(1) = nl
	sz(2) = nl
	
	stt = DftiCreateDescriptor( MD_Handle, DFTI_SINGLE, DFTI_COMPLEX, 2, sz)
	stt = DftiSetValue( MD_Handle, DFTI_PLACEMENT, DFTI_INPLACE)
	stt = DftiCommitDescriptor( MD_Handle)
	stt = DftiComputeForward( MD_Handle, sitecorr)
	stt = DftiFreeDescriptor( MD_Handle)
	
end subroutine SiteCorrelator

subroutine FourierForAllSites( nl, sits, zout)
	! Variables
	implicit NONE
	integer(4), intent(IN)			:: nl
	integer(4), intent(IN)			:: sits(nl,nl)
	complex, intent(OUT)			:: zout(nl*nl)
!---------------------------------------------------------------
	integer(4)					:: sts
	integer(4)					:: ix, iy
	complex, allocatable		:: ztmp(:)

	! Body of FourierForAllSites
	
	allocate( ztmp(nl*nl), stat = sts)
	zout = 0.0d0
	
	do ix = 1, nl
		do iy = 1, nl
			call SiteCorrelator( nl, ix, iy, sits, ztmp)
			zout = zout + ztmp
		end do
	end do
	
	zout = zout / ( nl*nl)
	
	deallocate( ztmp)

end subroutine FourierForAllSites

end module fouriers2d_m